library("mediation")
library('data.table')
library("dplyr")
library("dotwhisker")
library("ggpubr")
library("mma")

## ----------------------
## Main Result
## ----------------------

main.result <- function(data1, data2, plot.covs, title.l, title.r) {
  res1_pti_1 <- res2_pti_1 <- res3_pti_1 <- res4_pti_1 <- list()
  res1_pti_se_1 <- res2_pti_se_1 <- res3_pti_se_1 <- res4_pti_se_1 <- list()
  
  res1_pti_2 <- res2_pti_2 <- res3_pti_2 <- res4_pti_2 <- list()
  res1_pti_se_2 <- res2_pti_se_2 <- res3_pti_se_2 <- res4_pti_se_2 <- list()
  
  for(i in 1:3) {
    data.sub <- data1[data1$Group %in% c(0, i), ]    
    
    main1_1 <- lm(Military ~ as.factor(RussianGroup), data = data.sub)
    main2_1 <- lm(Military ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    main3_1 <- lm(Taiwan ~ as.factor(RussianGroup), data = data.sub)
    main4_1 <- lm(Taiwan ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    
    data.sub <- data2[data2$Group %in% c(0, i), ]    
    
    main1_2 <- lm(Military ~ as.factor(RussianGroup), data = data.sub)
    main2_2 <- lm(Military ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    main3_2 <- lm(Taiwan ~ as.factor(RussianGroup), data = data.sub)
    main4_2 <- lm(Taiwan ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    
    res1_pti_1[[i]] <- coefficients(main1_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_1)))]
    res1_pti_se_1[[i]] <- coef(summary(main1_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_1))), "Std. Error"]
    res3_pti_1[[i]] <- coefficients(main3_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_1)))]
    res3_pti_se_1[[i]] <- coef(summary(main3_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_1))), "Std. Error"]
    
    res2_pti_1[[i]] <- coefficients(main2_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_1)))]
    res2_pti_se_1[[i]] <- coef(summary(main2_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_1))), "Std. Error"]        
    res4_pti_1[[i]] <- coefficients(main4_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_1)))]
    res4_pti_se_1[[i]] <- coef(summary(main4_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_1))), "Std. Error"]    
    
    res1_pti_2[[i]] <- coefficients(main1_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_2)))]
    res1_pti_se_2[[i]] <- coef(summary(main1_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_2))), "Std. Error"]
    res3_pti_2[[i]] <- coefficients(main3_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_2)))]
    res3_pti_se_2[[i]] <- coef(summary(main3_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_2))), "Std. Error"]
    
    res2_pti_2[[i]] <- coefficients(main2_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_2)))]
    res2_pti_se_2[[i]] <- coef(summary(main2_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_2))), "Std. Error"]        
    res4_pti_2[[i]] <- coefficients(main4_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_2)))]
    res4_pti_se_2[[i]] <- coef(summary(main4_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_2))), "Std. Error"]    
  }
  
  if(plot.covs == TRUE) {
    m1_1 <- data.frame(cbind(do.call('rbind', res1_pti_1), do.call('rbind', res1_pti_se_1)))
    m3_1 <- data.frame(cbind(do.call('rbind', res3_pti_1), do.call('rbind', res3_pti_se_1)))
    
    m1_2 <- data.frame(cbind(do.call('rbind', res1_pti_2), do.call('rbind', res1_pti_se_2)))
    m3_2 <- data.frame(cbind(do.call('rbind', res3_pti_2), do.call('rbind', res3_pti_se_2)))
    
    m1_1$term <- paste0("t", 1:nrow(m1_1))
    m3_1$term <- paste0("t", 1:nrow(m3_1))
    
    m1_2$term <- paste0("t", 1:nrow(m1_2))
    m3_2$term <- paste0("t", 1:nrow(m3_2))
    
    colnames(m1_1) <- colnames(m3_1) <- c("estimate", "std.error", "term")
    colnames(m1_2) <- colnames(m3_2) <- c("estimate", "std.error", "term")
    
    m1_1$model <- 0
    m1_2$model <- 1
    
    m1 <- rbind(m1_2, m1_1)
    
    g1 <- {dwplot(m1, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for the Use of Force") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle(title.l) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23),
              legend.text = element_text(size = 23),
              legend.position = c(.43, 1.10),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))} 
    m3_1$model <- 0
    m3_2$model <- 1
    
    m3 <- rbind(m3_2, m3_1)
    
    g3 <- {dwplot(m3, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for Use of Force Against Taiwan") + ylab("Treatment Condition") + xlim(-0.3, 0.4) +
        ggtitle(title.r) +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) + 
        scale_colour_grey(start = 0.1, end = .7)} 
    
    ggarrange(g1, g3, ncol = 2, nrow = 1)
  } else {
    m2_1 <- data.frame(cbind(do.call('rbind', res2_pti_1), do.call('rbind', res2_pti_se_1)))
    m4_1 <- data.frame(cbind(do.call('rbind', res4_pti_1), do.call('rbind', res4_pti_se_1)))
    
    m2_2 <- data.frame(cbind(do.call('rbind', res2_pti_2), do.call('rbind', res2_pti_se_2)))
    m4_2 <- data.frame(cbind(do.call('rbind', res4_pti_2), do.call('rbind', res4_pti_se_2)))
    
    m2_1$term <- paste0("t", 1:nrow(m2_1))
    m4_1$term <- paste0("t", 1:nrow(m4_1))
    
    m2_2$term <- paste0("t", 1:nrow(m2_2))
    m4_2$term <- paste0("t", 1:nrow(m4_2))
    
    colnames(m2_1) <- colnames(m2_2) <- c("estimate", "std.error", "term")
    colnames(m4_1) <- colnames(m4_2) <- c("estimate", "std.error", "term")
    
    m2_1$model <- 0
    m2_2$model <- 1
    
    m2 <- rbind(m2_2, m2_1)
    
    g2 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for the Use of Force") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle(title.l) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23),
              legend.text = element_text(size = 23),
              legend.position = c(.43, 1.10),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))} 
    m4_1$model <- 0
    m4_2$model <- 1
    
    m4 <- rbind(m4_2, m4_1)
    
    g4 <- {dwplot(m4, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for Use of Force Against Taiwan") + ylab("Treatment Condition") + xlim(-0.3, 0.4) +
        ggtitle(title.r) +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) + 
        scale_colour_grey(start = 0.1, end = .7)} 
    
    ggarrange(g2, g4, ncol = 2, nrow = 1)
  }
}

## ----------------------
## Main Result: Different Baseline
## ----------------------

main.result.inva <- function(data1, data2, plot.covs,  title.l, title.r) {
  res1_pti_1 <- res2_pti_1 <- res3_pti_1 <- res4_pti_1 <- list()
  res1_pti_se_1 <- res2_pti_se_1 <- res3_pti_se_1 <- res4_pti_se_1 <- list()
  
  res1_pti_2 <- res2_pti_2 <- res3_pti_2 <- res4_pti_2 <- list()
  res1_pti_se_2 <- res2_pti_se_2 <- res3_pti_se_2 <- res4_pti_se_2 <- list()
  
  for(i in 1:2) {
    data.sub <- data1[!(data1$Group %in% 0), ]    
    data.sub <- data1[data1$Group %in% c(1, i + 1), ]    
    
    main1_1 <- lm(Military ~ as.factor(RussianGroup), data = data.sub)
    main2_1 <- lm(Military ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    main3_1 <- lm(Taiwan ~ as.factor(RussianGroup), data = data.sub)
    main4_1 <- lm(Taiwan ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    
    data.sub <- data2[!(data2$Group %in% 0), ]    
    data.sub <- data2[data2$Group %in% c(1, i + 1), ]  
    
    main1_2 <- lm(Military ~ as.factor(RussianGroup), data = data.sub)
    main2_2 <- lm(Military ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    main3_2 <- lm(Taiwan ~ as.factor(RussianGroup), data = data.sub)
    main4_2 <- lm(Taiwan ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    
    res1_pti_1[[i]] <- coefficients(main1_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_1)))]
    res1_pti_se_1[[i]] <- coef(summary(main1_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_1))), "Std. Error"]
    res3_pti_1[[i]] <- coefficients(main3_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_1)))]
    res3_pti_se_1[[i]] <- coef(summary(main3_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_1))), "Std. Error"]
    
    res2_pti_1[[i]] <- coefficients(main2_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_1)))]
    res2_pti_se_1[[i]] <- coef(summary(main2_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_1))), "Std. Error"]        
    res4_pti_1[[i]] <- coefficients(main4_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_1)))]
    res4_pti_se_1[[i]] <- coef(summary(main4_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_1))), "Std. Error"]    
    
    res1_pti_2[[i]] <- coefficients(main1_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_2)))]
    res1_pti_se_2[[i]] <- coef(summary(main1_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_2))), "Std. Error"]
    res3_pti_2[[i]] <- coefficients(main3_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_2)))]
    res3_pti_se_2[[i]] <- coef(summary(main3_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_2))), "Std. Error"]
    
    res2_pti_2[[i]] <- coefficients(main2_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_2)))]
    res2_pti_se_2[[i]] <- coef(summary(main2_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_2))), "Std. Error"]        
    res4_pti_2[[i]] <- coefficients(main4_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_2)))]
    res4_pti_se_2[[i]] <- coef(summary(main4_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_2))), "Std. Error"]    
  }
  
  if(plot.covs == TRUE) {
    m1_1 <- data.frame(cbind(do.call('rbind', res1_pti_1), do.call('rbind', res1_pti_se_1)))
    m3_1 <- data.frame(cbind(do.call('rbind', res3_pti_1), do.call('rbind', res3_pti_se_1)))
    
    m1_2 <- data.frame(cbind(do.call('rbind', res1_pti_2), do.call('rbind', res1_pti_se_2)))
    m3_2 <- data.frame(cbind(do.call('rbind', res3_pti_2), do.call('rbind', res3_pti_se_2)))
    
    m1_1$term <- paste0("t", 1:nrow(m1_1))
    m3_1$term <- paste0("t", 1:nrow(m3_1))
    
    m1_2$term <- paste0("t", 1:nrow(m1_2))
    m3_2$term <- paste0("t", 1:nrow(m3_2))
    
    colnames(m1_1) <- colnames(m3_1) <- c("estimate", "std.error", "term")
    colnames(m1_2) <- colnames(m3_2) <- c("estimate", "std.error", "term")
    
    m1_1$model <- 0
    m1_2$model <- 1
    
    m1 <- rbind(m1_2, m1_1)
    
    g1 <- {dwplot(m1, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Economic\n Measures", t2 = "Military\n Measures", t3 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for the Use of Force") + ylab("Treatment Condition") + xlim(-0.4, 0.4)+
        ggtitle(title.l) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23),
              legend.text = element_text(size = 23),
              legend.position = c(.43, 1.10),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))} 
    m3_1$model <- 0
    m3_2$model <- 1
    
    m3 <- rbind(m3_2, m3_1)
    
    g3 <- {dwplot(m3, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Economic\n Measures", t2 = "Military\n Measures", t3 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for Use of Force Against Taiwan") + ylab("Treatment Condition") + xlim(-0.4, 0.4) +
        ggtitle(title.r) +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) + 
        scale_colour_grey(start = 0.1, end = .7)} 
    
    ggarrange(g1, g3, ncol = 2, nrow = 1)
  } else {
    m2_1 <- data.frame(cbind(do.call('rbind', res2_pti_1), do.call('rbind', res2_pti_se_1)))
    m4_1 <- data.frame(cbind(do.call('rbind', res4_pti_1), do.call('rbind', res4_pti_se_1)))
    
    m2_2 <- data.frame(cbind(do.call('rbind', res2_pti_2), do.call('rbind', res2_pti_se_2)))
    m4_2 <- data.frame(cbind(do.call('rbind', res4_pti_2), do.call('rbind', res4_pti_se_2)))
    
    m2_1$term <- paste0("t", 1:nrow(m2_1))
    m4_1$term <- paste0("t", 1:nrow(m4_1))
    
    m2_2$term <- paste0("t", 1:nrow(m2_2))
    m4_2$term <- paste0("t", 1:nrow(m4_2))
    
    colnames(m2_1) <- colnames(m2_2) <- c("estimate", "std.error", "term")
    colnames(m4_1) <- colnames(m4_2) <- c("estimate", "std.error", "term")
    
    m2_1$model <- 0
    m2_2$model <- 1
    
    m2 <- rbind(m2_2, m2_1)
    
    g2 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Economic\n Measures", t2 = "Military\n Measures", t3 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for the Use of Force") + ylab("Treatment Condition") + xlim(-0.4, 0.4)+
        ggtitle(title.l) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23),
              legend.text = element_text(size = 23),
              legend.position = c(.43, 1.10),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))} 
    m4_1$model <- 0
    m4_2$model <- 1
    
    m4 <- rbind(m4_2, m4_1)
    
    g4 <- {dwplot(m4, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Economic\n Measures", t2 = "Military\n Measures", t3 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for Use of Force Against Taiwan") + ylab("Treatment Condition") + xlim(-0.4, 0.4) +
        ggtitle(title.r) +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) + 
        scale_colour_grey(start = 0.1, end = .7)} 
    
    ggarrange(g2, g4, ncol = 2, nrow = 1)
  }
}

## ----------------------
## Main Result: New Taiwan Question
## ----------------------

main.result.taiwan <- function(data, plot.covs, paper,  title.l, title.r) {
  res1_pti <- res2_pti <- res3_pti <- res4_pti <- res5_pti <- res6_pti <- res7_pti <- res8_pti <- res9_pti <- res10_pti <-list()
  res1_pti_se <- res2_pti_se <- res3_pti_se <- res4_pti_se <- res5_pti_se <- res6_pti_se <- res7_pti_se <- res8_pti_se <- res9_pti_se <- res10_pti_se <- list()
  
  for(i in 1:3) {
    data.sub <- data[data$Group %in% c(0, i), ]    
    
    main1 <- lm(Taiwan_Invasion ~ as.factor(RussianGroup), data = data.sub)
    main2 <- lm(Taiwan_Invasion ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                data = data.sub)
    main3 <- lm(Taiwan_Coersion ~ as.factor(RussianGroup), data = data.sub)
    main4 <- lm(Taiwan_Coersion ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                data = data.sub)
    main5 <- lm(Taiwan_Sanction ~ as.factor(RussianGroup), data = data.sub)
    main6 <- lm(Taiwan_Sanction ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                data = data.sub)
    main7 <- lm(Taiwan_StatusQuo ~ as.factor(RussianGroup), data = data.sub)
    main8 <- lm(Taiwan_StatusQuo ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                data = data.sub)
    main9 <- lm(Taiwan_Independent ~ as.factor(RussianGroup), data = data.sub)
    main10 <- lm(Taiwan_Independent ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                 + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                 + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                 data = data.sub)
    
    res1_pti[[i]] <- coefficients(main1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1)))]
    res1_pti_se[[i]] <- coef(summary(main1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1))), "Std. Error"]
    res3_pti[[i]] <- coefficients(main3)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3)))]
    res3_pti_se[[i]] <- coef(summary(main3))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3))), "Std. Error"]
    res5_pti[[i]] <- coefficients(main5)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main5)))]
    res5_pti_se[[i]] <- coef(summary(main5))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4))), "Std. Error"]
    res7_pti[[i]] <- coefficients(main7)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main7)))]
    res7_pti_se[[i]] <- coef(summary(main7))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main7))), "Std. Error"]
    res9_pti[[i]] <- coefficients(main9)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main9)))]
    res9_pti_se[[i]] <- coef(summary(main9))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main9))), "Std. Error"]
    
    res2_pti[[i]] <- coefficients(main2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2)))]
    res2_pti_se[[i]] <- coef(summary(main2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2))), "Std. Error"]        
    res4_pti[[i]] <- coefficients(main4)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4)))]
    res4_pti_se[[i]] <- coef(summary(main4))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4))), "Std. Error"]    
    res6_pti[[i]] <- coefficients(main6)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main6)))]
    res6_pti_se[[i]] <- coef(summary(main6))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main6))), "Std. Error"]    
    res8_pti[[i]] <- coefficients(main8)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main8)))]
    res8_pti_se[[i]] <- coef(summary(main8))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main8))), "Std. Error"]    
    res10_pti[[i]] <- coefficients(main10)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main10)))]
    res10_pti_se[[i]] <- coef(summary(main10))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main10))), "Std. Error"]    
    
  }
  
  if(plot.covs == TRUE) {
    m1 <- data.frame(cbind(do.call('rbind', res1_pti), do.call('rbind', res1_pti_se)))
    m3 <- data.frame(cbind(do.call('rbind', res3_pti), do.call('rbind', res3_pti_se)))
    m5 <- data.frame(cbind(do.call('rbind', res5_pti), do.call('rbind', res5_pti_se)))
    m7 <- data.frame(cbind(do.call('rbind', res7_pti), do.call('rbind', res7_pti_se)))
    m9 <- data.frame(cbind(do.call('rbind', res9_pti), do.call('rbind', res9_pti_se)))
    
    
    m1$term <- paste0("t", 1:nrow(m1))
    m3$term <- paste0("t", 1:nrow(m3))
    m5$term <- paste0("t", 1:nrow(m5))
    m7$term <- paste0("t", 1:nrow(m7))
    m9$term <- paste0("t", 1:nrow(m9))
    
    colnames(m1) <- colnames(m3) <- colnames(m5) <- colnames(m7) <- colnames(m9) <- c("estimate", "std.error", "term")
    
    g1 <- {dwplot(m1, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Unification via War") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle(title.l) +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) +
        scale_colour_grey(start = 0.1, end = .7)} 
    
    g3 <- {dwplot(m3, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Unification via Military Coercion") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle(title.r) +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) +
        scale_colour_grey(start = 0.1, end = .7)} 
    
    g5 <- {dwplot(m5, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Unification via Economic Sanctions") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) +
        scale_colour_grey(start = 0.1, end = .7)} 
    
    g7 <- {dwplot(m7, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Unification via Status Quo") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23))+
        scale_colour_grey(start = 0.1, end = .7)} 
    
    g9 <- {dwplot(m9, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Taiwan's Independence") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23))+
        scale_colour_grey(start = 0.1, end = .7)} 
    
    if(paper == F) {
      ggarrange(g1, g3, g5, g7, g9, ncol = 3, nrow = 2)
    } else {
      ggarrange(g1, g3, ncol = 2, nrow = 1)
    }
    
  } else {
    m2 <- data.frame(cbind(do.call('rbind', res2_pti), do.call('rbind', res2_pti_se)))
    m4 <- data.frame(cbind(do.call('rbind', res4_pti), do.call('rbind', res4_pti_se)))
    m6 <- data.frame(cbind(do.call('rbind', res6_pti), do.call('rbind', res6_pti_se)))
    m8 <- data.frame(cbind(do.call('rbind', res8_pti), do.call('rbind', res8_pti_se)))
    m10 <- data.frame(cbind(do.call('rbind', res10_pti), do.call('rbind', res10_pti_se)))
    
    m2$term <- paste0("t", 1:nrow(m2))
    m4$term <- paste0("t", 1:nrow(m4))
    m6$term <- paste0("t", 1:nrow(m6))
    m8$term <- paste0("t", 1:nrow(m8))
    m10$term <- paste0("t", 1:nrow(m10))
    
    colnames(m2) <- colnames(m4) <- colnames(m6) <- colnames(m8) <- colnames(m10) <- c("estimate", "std.error", "term")
    
    g2 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Unification via War") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle(title.l) +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23))} 
    
    g4 <- {dwplot(m4, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Taiwan Coercion") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle(title.r) +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23))} 
    
    g6 <- {dwplot(m6, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Unification via Military Coercion") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23))} 
    
    g8 <- {dwplot(m8, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Taiwan Status Quo") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23))} 
    
    g10 <- {dwplot(m10, 
                   dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                   vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Taiwan Independence") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23))} 
    if(paper == F) {
      ggarrange(g2, g4, g6, g8, g10, ncol = 3, nrow = 2)
    } else {
      ggarrange(g2, g4, ncol = 2, nrow = 1)
    }
  }
}

## ----------------------
## Mediation with Multiple Mediators
## ----------------------

mediation.multAll <- function(data1, data2, treatment.value, control.value, outcome, alt = FALSE, title) {
  
  alp <- .20
  
  data.sub <- data1[data1$RussianGroup %in% c(control.value, treatment.value), ]
  x <- data.sub[, c("Econ_Cost", "Military_Cost", "Success", "Threat", "Morality", "US", "Peace"), with = F]
  colnames(x) <- c("Economic_cost", "Military_cost", "Success", "Threat", "Morality", "US", "Peaceful")
  pred <- data.sub[, "Group"]
  y <- as.matrix(data.sub[, outcome, with = F])
  data.1 <- data.org(x, y, mediator = 1:7,
                     pred = pred,
                     alpha=alp, alpha2=alp)
  
  med.1 <- boot.med(data = data.1, n2 = 500)
  out1 <- summary(med.1)
  xlab <- "Indirect Effect"
  
  data.sub <- data2[data2$RussianGroup %in% c(control.value, treatment.value), ]
  
  x <- data.sub[, c("Econ_Cost", "Military_Cost", "Success", "Threat", "Morality", "US", "Peace", "Legality", "US_Threat", "Peaceful_Image"), with = F]
  colnames(x) <- c("Economic_cost", "Military_cost", "Success", "Threat", "Morality", "US", "Peaceful", "Legality", "US_Threat", "Peaceful_Image")
  pred <- data.sub[, "Group"]
  y <- as.matrix(data.sub[, outcome, with = F])
  data.2 <- data.org(x, y, mediator = 1:10,
                     pred = pred,
                     alpha=alp, alpha2=alp)
  
  med.2 <- boot.med(data = data.2, n2 = 500)
  out2 <- summary(med.2)
  
  
  if(treatment.value == "Invasion" & alt == FALSE) {
    treat.name <- "A. Invasion"
  }
  
  if(treatment.value == "Invasion" & alt == TRUE & outcome == "Military") {
    treat.name <- "Support for\nUse of Force"
    xlab <- "Indirect Effect of Invasion"
  }
  
  if(treatment.value == "Invasion" & alt == TRUE & outcome == "Taiwan") {
    treat.name <- "Support for Use of\nForce Against Taiwan"
    xlab <- "Indirect Effect of Invasion"
  }
  
  if(treatment.value == "Sanction") {
    treat.name <- "B. Economic Sanction"
  }
  
  if(treatment.value == "Military") {
    treat.name <- "C. Military Action"
  }
  
  if(treatment.value == "LackMilitary") {
    treat.name <- "D. Lack of Military Action"
  }
  
  res2_1 <- out1$bin.result$results$indirect.effect$Group["est", ]
  res2_1 <- res2_1[names(res2_1) %in% c("y1.all", "y1.Success", "y1.Morality", "y1.Peaceful")]
  
  res2_2 <- out2$bin.result$results$indirect.effect$Group["est", ]
  res2_2 <- res2_2[names(res2_2) %in% c("y1.all", "y1.Success", "y1.Morality", "y1.Peaceful")]
  
  res2_se_1 <- out1$bin.result$results$indirect.effect$Group["sd", ]
  res2_se_1 <- res2_se_1[names(res2_se_1) %in% c("y1.all", "y1.Success", "y1.Morality", "y1.Peaceful")]
  
  res2_se_2 <- out2$bin.result$results$indirect.effect$Group["sd", ]
  res2_se_2 <- res2_se_2[names(res2_se_2) %in% c("y1.all", "y1.Success", "y1.Morality", "y1.Peaceful")]
  
  m2_1 <- data.frame(cbind(res2_1, res2_se_1))
  order <- order(m2_1$res2_1, decreasing = T)
  m2_1 <- m2_1[order, ]
  
  m2_2 <- data.frame(cbind(res2_2, res2_se_2))
  m2_2 <- m2_2[order, ]
  
  m2_1$term <- paste0("t", 1:nrow(m2_1))
  m2_2$term <- paste0("t", 1:nrow(m2_2))
  
  colnames(m2_1) <- c("estimate", "std.error", "term")
  rownames(m2_1) <- gsub("y1.", "", rownames(m2_1))
  
  colnames(m2_2) <- c("estimate", "std.error", "term")
  rownames(m2_2) <- gsub("y1.", "", rownames(m2_2))
  
  m2_1$model <- 0
  m2_2$model <- 1
  
  m2 <- rbind(m2_2, m2_1)
  
  if(outcome == "Military") {
    g1 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "All Mediators", 
                             t2 = rownames(m2)[2], 
                             t3 = rownames(m2)[3], 
                             t4 = rownames(m2)[4])) +
        theme_classic() + xlab(xlab) + ylab("Mediator") + xlim(-0.15, 0.30)+
        ggtitle(paste0(title, treat.name)) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 18),
              axis.text.y  = element_text(size = 18),
              text = element_text(size = 18),
              legend.text = element_text(size = 18),
              legend.position = c(1.0, 0.22),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))      } } else {
                            g1 <- {dwplot(m2, 
                                          dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                                          vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
                                relabel_predictors(c(t1 = "All Mediators", 
                                                     t2 = rownames(m2)[2], 
                                                     t3 = rownames(m2)[3], 
                                                     t4 = rownames(m2)[4])) +
                                theme_classic() + xlab(xlab) + ylab("Mediator") + xlim(-0.15, 0.30)+
                                ggtitle(paste0(title, treat.name)) +
                                theme(plot.title = element_text(face="bold"), 
                                      axis.text.x  = element_text(size = 18),
                                      axis.text.y  = element_text(size = 18),
                                      text = element_text(size = 18), legend.position = "none") + 
                                scale_colour_grey(start = 0.1, end = .7)      }                          
                            
                          }
  g1
}

## ----------------------
## Main Result Appendix
## ----------------------

main.result.app <- function(data1, data2, plot.covs) {
  res1_pti_1 <- res2_pti_1 <- res3_pti_1 <- res4_pti_1 <- list()
  res1_pti_se_1 <- res2_pti_se_1 <- res3_pti_se_1 <- res4_pti_se_1 <- list()
  
  res1_pti_2 <- res2_pti_2 <- res3_pti_2 <- res4_pti_2 <- list()
  res1_pti_se_2 <- res2_pti_se_2 <- res3_pti_se_2 <- res4_pti_se_2 <- list()
  
  for(i in 1:4) {
    data.sub <- data1[data1$Group %in% c(0, i), ]    
    
    main1_1 <- lm(Military ~ as.factor(RussianGroup), data = data.sub)
    main2_1 <- lm(Military ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    main3_1 <- lm(Taiwan ~ as.factor(RussianGroup), data = data.sub)
    main4_1 <- lm(Taiwan ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    
    res1_pti_1[[i]] <- coefficients(main1_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_1)))]
    res1_pti_se_1[[i]] <- coef(summary(main1_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_1))), "Std. Error"]
    res3_pti_1[[i]] <- coefficients(main3_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_1)))]
    res3_pti_se_1[[i]] <- coef(summary(main3_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_1))), "Std. Error"]
    
    res2_pti_1[[i]] <- coefficients(main2_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_1)))]
    res2_pti_se_1[[i]] <- coef(summary(main2_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_1))), "Std. Error"]        
    res4_pti_1[[i]] <- coefficients(main4_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_1)))]
    res4_pti_se_1[[i]] <- coef(summary(main4_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_1))), "Std. Error"]    
    
    if(i != 4) {
      data.sub <- data2[data2$Group %in% c(0, i), ]    
      
      main1_2 <- lm(Military ~ as.factor(RussianGroup), data = data.sub)
      main2_2 <- lm(Military ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                    + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                    + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                    data = data.sub)
      main3_2 <- lm(Taiwan ~ as.factor(RussianGroup), data = data.sub)
      main4_2 <- lm(Taiwan ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                    + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                    + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                    data = data.sub)
      
      res1_pti_2[[i]] <- coefficients(main1_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_2)))]
      res1_pti_se_2[[i]] <- coef(summary(main1_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_2))), "Std. Error"]
      res3_pti_2[[i]] <- coefficients(main3_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_2)))]
      res3_pti_se_2[[i]] <- coef(summary(main3_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_2))), "Std. Error"]
      
      res2_pti_2[[i]] <- coefficients(main2_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_2)))]
      res2_pti_se_2[[i]] <- coef(summary(main2_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_2))), "Std. Error"]        
      res4_pti_2[[i]] <- coefficients(main4_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_2)))]
      res4_pti_se_2[[i]] <- coef(summary(main4_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_2))), "Std. Error"]    }
    else {
      res1_pti_2[[i]] <- NA
      res1_pti_se_2[[i]] <- NA
      res3_pti_2[[i]] <- NA
      res3_pti_se_2[[i]] <- NA
      
      res2_pti_2[[i]] <- NA
      res2_pti_se_2[[i]] <- NA
      res4_pti_2[[i]] <- NA
      res4_pti_se_2[[i]] <- NA
    }
  }
  
  if(plot.covs == TRUE) {
    m1_1 <- data.frame(cbind(do.call('rbind', res1_pti_1), do.call('rbind', res1_pti_se_1)))
    m3_1 <- data.frame(cbind(do.call('rbind', res3_pti_1), do.call('rbind', res3_pti_se_1)))
    
    m1_2 <- data.frame(cbind(do.call('rbind', res1_pti_2), do.call('rbind', res1_pti_se_2)))
    m3_2 <- data.frame(cbind(do.call('rbind', res3_pti_2), do.call('rbind', res3_pti_se_2)))
    
    m1_1$term <- paste0("t", 1:nrow(m1_1))
    m3_1$term <- paste0("t", 1:nrow(m3_1))
    
    m1_2$term <- paste0("t", 1:nrow(m1_2))
    m3_2$term <- paste0("t", 1:nrow(m3_2))
    
    colnames(m1_1) <- colnames(m3_1) <- c("estimate", "std.error", "term")
    colnames(m1_2) <- colnames(m3_2) <- c("estimate", "std.error", "term")
    
    m1_1$model <- 0
    m1_2$model <- 1
    
    m1 <- rbind(m1_2, m1_1)
    
    g1 <- {dwplot(m1, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for the Use of Force") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23),
              legend.text = element_text(size = 23),
              legend.position = c(1.00, 0.28),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))} 
    m3_1$model <- 0
    m3_2$model <- 1
    
    m3 <- rbind(m3_2, m3_1)
    
    g3 <- {dwplot(m3, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for Use of Force Against Taiwan") + ylab("Treatment Condition") + xlim(-0.3, 0.4) +
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) + 
        scale_colour_grey(start = 0.1, end = .7)} 
    
    ggarrange(g1, g3, ncol = 2, nrow = 1)
  } else {
    m2_1 <- data.frame(cbind(do.call('rbind', res2_pti_1), do.call('rbind', res2_pti_se_1)))
    m4_1 <- data.frame(cbind(do.call('rbind', res4_pti_1), do.call('rbind', res4_pti_se_1)))
    
    m2_2 <- data.frame(cbind(do.call('rbind', res2_pti_2), do.call('rbind', res2_pti_se_2)))
    m4_2 <- data.frame(cbind(do.call('rbind', res4_pti_2), do.call('rbind', res4_pti_se_2)))
    
    m2_1$term <- paste0("t", 1:nrow(m2_1))
    m4_1$term <- paste0("t", 1:nrow(m4_1))
    
    m2_2$term <- paste0("t", 1:nrow(m2_2))
    m4_2$term <- paste0("t", 1:nrow(m4_2))
    
    colnames(m2_1) <- colnames(m2_2) <- c("estimate", "std.error", "term")
    colnames(m4_1) <- colnames(m4_2) <- c("estimate", "std.error", "term")
    
    m2_1$model <- 0
    m2_2$model <- 1
    
    m2 <- rbind(m2_2, m2_1)
    
    g2 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for the Use of Force") + ylab("Treatment Condition") + xlim(-0.3, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23),
              legend.text = element_text(size = 23),
              legend.position = c(1.00, 0.28),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))} 
    m4_1$model <- 0
    m4_2$model <- 1
    
    m4 <- rbind(m4_2, m4_1)
    
    g4 <- {dwplot(m4, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Invasion", t2 = "Economic\n Measures", t3 = "Military\n Measures", t4 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for Use of Force Against Taiwan") + ylab("Treatment Condition") + xlim(-0.3, 0.4) +
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) + 
        scale_colour_grey(start = 0.1, end = .7)} 
    
    ggarrange(g2, g4, ncol = 2, nrow = 1)
  }
}

## ----------------------
## Main Result Different Baseline Appendix
## ----------------------

main.result.inva.app <- function(data1, data2, plot.covs) {
  res1_pti_1 <- res2_pti_1 <- res3_pti_1 <- res4_pti_1 <- list()
  res1_pti_se_1 <- res2_pti_se_1 <- res3_pti_se_1 <- res4_pti_se_1 <- list()
  
  res1_pti_2 <- res2_pti_2 <- res3_pti_2 <- res4_pti_2 <- list()
  res1_pti_se_2 <- res2_pti_se_2 <- res3_pti_se_2 <- res4_pti_se_2 <- list()
  
  for(i in 1:3) {
    data.sub <- data1[!(data1$Group %in% 0), ]    
    data.sub <- data1[data1$Group %in% c(1, i + 1), ]    
    
    main1_1 <- lm(Military ~ as.factor(RussianGroup), data = data.sub)
    main2_1 <- lm(Military ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    main3_1 <- lm(Taiwan ~ as.factor(RussianGroup), data = data.sub)
    main4_1 <- lm(Taiwan ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                  + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                  + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                  data = data.sub)
    
    res1_pti_1[[i]] <- coefficients(main1_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_1)))]
    res1_pti_se_1[[i]] <- coef(summary(main1_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_1))), "Std. Error"]
    res3_pti_1[[i]] <- coefficients(main3_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_1)))]
    res3_pti_se_1[[i]] <- coef(summary(main3_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_1))), "Std. Error"]
    
    res2_pti_1[[i]] <- coefficients(main2_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_1)))]
    res2_pti_se_1[[i]] <- coef(summary(main2_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_1))), "Std. Error"]        
    res4_pti_1[[i]] <- coefficients(main4_1)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_1)))]
    res4_pti_se_1[[i]] <- coef(summary(main4_1))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_1))), "Std. Error"]    
    
    if(i != 3) {
      data.sub <- data2[!(data2$Group %in% 0), ]    
      data.sub <- data2[data2$Group %in% c(1, i + 1), ]  
      
      main1_2 <- lm(Military ~ as.factor(RussianGroup), data = data.sub)
      main2_2 <- lm(Military ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                    + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                    + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                    data = data.sub)
      main3_2 <- lm(Taiwan ~ as.factor(RussianGroup), data = data.sub)
      main4_2 <- lm(Taiwan ~ as.factor(RussianGroup) + as.factor(Female) + as.factor(Age_Group)
                    + as.factor(Education) + as.factor(PartyMember) + as.factor(Pol_Interest)
                    + as.factor(Ideology) + as.factor(Nationalism) + as.factor(Social_Media) + as.factor(Foreign),
                    data = data.sub)
      
      res1_pti_2[[i]] <- coefficients(main1_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_2)))]
      res1_pti_se_2[[i]] <- coef(summary(main1_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main1_2))), "Std. Error"]
      res3_pti_2[[i]] <- coefficients(main3_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_2)))]
      res3_pti_se_2[[i]] <- coef(summary(main3_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main3_2))), "Std. Error"]
      
      res2_pti_2[[i]] <- coefficients(main2_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_2)))]
      res2_pti_se_2[[i]] <- coef(summary(main2_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main2_2))), "Std. Error"]        
      res4_pti_2[[i]] <- coefficients(main4_2)[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_2)))]
      res4_pti_se_2[[i]] <- coef(summary(main4_2))[grep("as.factor\\(RussianGroup\\)", names(coefficients(main4_2))), "Std. Error"]
    } else {
      res1_pti_2[[i]] <- NA
      res1_pti_se_2[[i]] <- NA
      res3_pti_2[[i]] <- NA
      res3_pti_se_2[[i]] <- NA
      
      res2_pti_2[[i]] <- NA
      res2_pti_se_2[[i]] <- NA
      res4_pti_2[[i]] <- NA
      res4_pti_se_2[[i]] <- NA
    }
  }
  
  if(plot.covs == TRUE) {
    m1_1 <- data.frame(cbind(do.call('rbind', res1_pti_1), do.call('rbind', res1_pti_se_1)))
    m3_1 <- data.frame(cbind(do.call('rbind', res3_pti_1), do.call('rbind', res3_pti_se_1)))
    
    m1_2 <- data.frame(cbind(do.call('rbind', res1_pti_2), do.call('rbind', res1_pti_se_2)))
    m3_2 <- data.frame(cbind(do.call('rbind', res3_pti_2), do.call('rbind', res3_pti_se_2)))
    
    m1_1$term <- paste0("t", 1:nrow(m1_1))
    m3_1$term <- paste0("t", 1:nrow(m3_1))
    
    m1_2$term <- paste0("t", 1:nrow(m1_2))
    m3_2$term <- paste0("t", 1:nrow(m3_2))
    
    colnames(m1_1) <- colnames(m3_1) <- c("estimate", "std.error", "term")
    colnames(m1_2) <- colnames(m3_2) <- c("estimate", "std.error", "term")
    
    m1_1$model <- 0
    m1_2$model <- 1
    
    m1 <- rbind(m1_2, m1_1)
    
    g1 <- {dwplot(m1, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Economic\n Measures", t2 = "Military\n Measures", t3 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for the Use of Force") + ylab("Treatment Condition") + xlim(-0.4, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23),
              legend.text = element_text(size = 23),
              legend.position = c(1.00, 0.28),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))} 
    m3_1$model <- 0
    m3_2$model <- 1
    
    m3 <- rbind(m3_2, m3_1)
    
    g3 <- {dwplot(m3, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Economic\n Measures", t2 = "Military\n Measures", t3 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for Use of Force Against Taiwan") + ylab("Treatment Condition") + xlim(-0.4, 0.4) +
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) + 
        scale_colour_grey(start = 0.1, end = .7)} 
    
    ggarrange(g1, g3, ncol = 2, nrow = 1)
  } else {
    m2_1 <- data.frame(cbind(do.call('rbind', res2_pti_1), do.call('rbind', res2_pti_se_1)))
    m4_1 <- data.frame(cbind(do.call('rbind', res4_pti_1), do.call('rbind', res4_pti_se_1)))
    
    m2_2 <- data.frame(cbind(do.call('rbind', res2_pti_2), do.call('rbind', res2_pti_se_2)))
    m4_2 <- data.frame(cbind(do.call('rbind', res4_pti_2), do.call('rbind', res4_pti_se_2)))
    
    m2_1$term <- paste0("t", 1:nrow(m2_1))
    m4_1$term <- paste0("t", 1:nrow(m4_1))
    
    m2_2$term <- paste0("t", 1:nrow(m2_2))
    m4_2$term <- paste0("t", 1:nrow(m4_2))
    
    colnames(m2_1) <- colnames(m2_2) <- c("estimate", "std.error", "term")
    colnames(m4_1) <- colnames(m4_2) <- c("estimate", "std.error", "term")
    
    m2_1$model <- 0
    m2_2$model <- 1
    
    m2 <- rbind(m2_2, m2_1)
    
    g2 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Economic\n Measures", t2 = "Military\n Measures", t3 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for the Use of Force") + ylab("Treatment Condition") + xlim(-0.4, 0.4)+
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23),
              legend.text = element_text(size = 23),
              legend.position = c(1.00, 0.28),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))} 
    m4_1$model <- 0
    m4_2$model <- 1
    
    m4 <- rbind(m4_2, m4_1)
    
    g4 <- {dwplot(m4, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Economic\n Measures", t2 = "Military\n Measures", t3 = "Lack of Military\n Measures")) +
        theme_classic() + xlab("The Effect of [Treatment Condition] on\n Support for Use of Force Against Taiwan") + ylab("Treatment Condition") + xlim(-0.4, 0.4) +
        ggtitle("") +
        theme(plot.title = element_text(face="bold"), legend.position = "none",
              axis.text.x  = element_text(size = 23),
              axis.text.y  = element_text(size = 23),
              text = element_text(size = 23)) + 
        scale_colour_grey(start = 0.1, end = .7)} 
    
    ggarrange(g2, g4, ncol = 2, nrow = 1)
  }
}


## ----------------------
## Mediation Decomposition
## ----------------------
mediation.mult.decomp <- function(data1, data2, treatment.value, control.value, outcome, alt = FALSE, legend) {
  
  alp <- 1
  data.sub <- data1[data1$RussianGroup %in% c(control.value, treatment.value), ]
  x <- data.sub[, c("Econ_Cost", "Military_Cost", "Success", "Threat", "Morality", "US", "Peace"), with = F]
  colnames(x) <- c("Economic_cost", "Military_cost", "Success", "Threat", "Morality", "US_Influence", "Peaceful")
  pred <- data.sub[, "Group"]
  y <- as.matrix(data.sub[, outcome, with = F])
  data.1 <- data.org(x, y, mediator = 1:7,
                     pred = pred,
                     alpha=alp, alpha2=alp)
  
  med.1 <- boot.med(data = data.1, n2 = 500)
  out1 <- summary(med.1)
  xlab <- "Indirect Effect"
  
  if(treatment.value != "LackMilitary") {
    data.sub <- data2[data2$RussianGroup %in% c(control.value, treatment.value), ]
    x <- data.sub[, c("Econ_Cost", "Military_Cost", "Success", "Threat", "Morality", "US", "Peace", "Legality", "US_Threat", "Peaceful_Image"), with = F]
    colnames(x) <- c("Economic_cost", "Military_cost", "Success", "Threat", "Morality", "US_Influence", "Peaceful", "Legality", "US_Threat", "Peaceful_Image")
    pred <- data.sub[, "Group"]
    y <- as.matrix(data.sub[, outcome, with = F])
    data.2 <- data.org(x, y, mediator = 1:10,
                       pred = pred,
                       alpha=alp, alpha2=alp)
    
    med.2 <- boot.med(data = data.2, n2 = 500)
    out2 <- summary(med.2)
  }
  
  if(outcome == "Military") {
    main1_1 <- lm(Military ~ as.factor(RussianGroup), data = data1)
    main1_2 <- lm(Military ~ as.factor(RussianGroup), data = data2)
  } else {
    main1_1 <- lm(Taiwan ~ as.factor(RussianGroup), data = data1)
    main1_2 <- lm(Taiwan ~ as.factor(RussianGroup), data = data2)
  }
  
  if(treatment.value == "Invasion" & alt == FALSE) {
    treat.name <- "A. Invasion"
  }
  
  if(treatment.value == "Invasion" & alt == TRUE & outcome == "Military") {
    treat.name <- ""
  }
  
  if(treatment.value == "Invasion" & alt == TRUE & outcome == "Taiwan") {
    treat.name <- ""
  }
  
  if(treatment.value == "Sanction") {
    treat.name <- "B. Economic Measures"
  }
  
  if(treatment.value == "Military") {
    treat.name <- "C. Military Measures"
  }
  
  if(treatment.value == "LackMilitary") {
    treat.name <- "D. Lack of Military Measures"
  }
  
  if(outcome == "Military") {
    low <- -0.15; up <- 0.30
  } else {
    low <- -0.15; up <- 0.35
  }
  
  # res2_1 <- round(out1$bin.result$results$indirect.effect$Group["est", ], 2)
  # res2_se_1 <- round(out1$bin.result$results$indirect.effect$Group["sd", ], 2)
  # 
  # res2_d <- round(out$bin.result$results$direct.effect["est", ], 2)
  # res2_se_d <- round(out$bin.result$results$direct.effect["sd", ], 2)
  # 
  # res2_t <- round(out$bin.result$results$total.effect["est", ], 2)
  # res2_se_t <- round(out$bin.result$results$total.effect["sd", ], 2)
  # 
  
  res2i_1 <- out1$bin.result$results$indirect.effect$Group["est", "y1.all"]
  res2i_se_1 <- out1$bin.result$results$indirect.effect$Group["sd", "y1.all"]
  res2d_1 <- out1$bin.result$results$direct.effect["est",]
  res2d_se_1 <- out1$bin.result$results$direct.effect["sd",]
  
  if(treatment.value != "LackMilitary") {
    res2i_2 <- out2$bin.result$results$indirect.effect$Group["est", "y1.all"]
    res2i_se_2 <- out2$bin.result$results$indirect.effect$Group["sd", "y1.all"]
    res2d_2 <- out2$bin.result$results$direct.effect["est",]
    res2d_se_2 <- out2$bin.result$results$direct.effect["sd",]
  }
  
  m2_1 <- data.frame(rbind(cbind(res2i_1, res2i_se_1), cbind(res2d_1, res2d_se_1)))
  
  if(treatment.value == "LackMilitary") {
    res2d_2 <- res2i_2 <- NA
    res2d_se_2 <- res2i_se_2 <- NA
    m2_2 <- data.frame(rbind(cbind(res2i_2, res2i_se_2), cbind(res2d_2, res2d_se_2)))
  } else {
    m2_2 <- data.frame(rbind(cbind(res2i_2, res2i_se_2), cbind(res2d_2, res2d_se_2)))
  }
  
  m2_1$term <- paste0("t", 1:nrow(m2_1))
  m2_2$term <- paste0("t", 1:nrow(m2_2))
  
  colnames(m2_1) <- c("estimate", "std.error", "term")
  colnames(m2_2) <- c("estimate", "std.error", "term")
  
  m2_1$model <- 0
  m2_2$model <- 1
  
  m2 <- rbind(m2_2, m2_1)
  #  rownames(m2) <- c("Direct Effect", "Indirect Effect\n (all mediators)", "Direct Effect", "Indirect Effect\n (all mediators)")
  
  if(outcome == "Military" & alt == FALSE) {
    xl <- "Effect of Treatment on\n Support for the Use of Force"
  } 
  
  if(outcome == "Military" & alt == TRUE) {
    xl <- "Effect of Invasion Treat. on\n Support for the Use of Force"
  } 
  
  if(outcome == "Taiwan" & alt == FALSE) {
    xl <- "Effect of Treatment on\n Support for the Use of Force Against Taiwan"
  } 
  
  if(outcome == "Taiwan" & alt == TRUE) {
    xl <- "Effect of Invasion Treat. on\n Support for the Use of Force Against Taiwan"
  } 
  
  if(legend == 0) {
    g1 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Direct Effect", 
                             t2 = "Indirect Effect\n (all mediators)")) +
        theme_classic() + xlab(xl) + ylab("") + xlim(c(low, up))+
        ggtitle(paste0("", treat.name)) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 18),
              axis.text.y  = element_text(size = 18),
              text = element_text(size = 18), legend.position = "none") + 
        scale_colour_grey(start = 0.1, end = .7)    } 
  } else {
    g1 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "Direct Effect", 
                             t2 = "Indirect Effect\n (all mediators)")) +
        theme_classic() + xlab(xl) + ylab("") + xlim(c(low, up))+
        ggtitle(paste0("", treat.name)) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 18),
              axis.text.y  = element_text(size = 18),
              text = element_text(size = 18),
              legend.text = element_text(size = 18),
              legend.position = c(1.0, 0.22),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))    } 
  }
  g1
}

## ----------------------
## Mediation with Multiple Mediators (all)
## ----------------------

mediation.multAll2 <- function(data1, data2, treatment.value, control.value, outcome, alt = FALSE, legend = 1) {
  
  alp <- 1
  
  data.sub <- data1[data1$RussianGroup %in% c(control.value, treatment.value), ]
  x <- data.sub[, c("Econ_Cost", "Military_Cost", "Success", "Threat", "Morality", "US", "Peace"), with = F]
  colnames(x) <- c("Economic_cost", "Military_cost", "Success", "Threat", "Morality", "US_Influence", "Peaceful")
  pred <- data.sub[, "Group"]
  y <- as.matrix(data.sub[, outcome, with = F])
  data.1 <- data.org(x, y, mediator = 1:7,
                     pred = pred,
                     alpha=alp, alpha2=alp)
  
  med.1 <- boot.med(data = data.1, n2 = 500)
  out1 <- summary(med.1)
  xlab <- "Indirect Effect"
  
  if(treatment.value != "LackMilitary") {
    data.sub <- data2[data2$RussianGroup %in% c(control.value, treatment.value), ]
    x <- data.sub[, c("Econ_Cost", "Military_Cost", "Success", "Threat", "Morality", "US", "Peace", "Legality", "US_Threat", "Peaceful_Image"), with = F]
    colnames(x) <- c("Economic_cost", "Military_cost", "Success", "Threat", "Morality", "US_Influence", "Peaceful", "Legality", "US_Threat", "Peaceful_Image")
    pred <- data.sub[, "Group"]
    y <- as.matrix(data.sub[, outcome, with = F])
    data.2 <- data.org(x, y, mediator = 1:10,
                       pred = pred,
                       alpha=alp, alpha2=alp)
    
    med.2 <- boot.med(data = data.2, n2 = 500)
    out2 <- summary(med.2)
  }
  
  if(treatment.value == "Invasion" & alt == FALSE) {
    treat.name <- "A. Invasion"
  }
  
  if(treatment.value == "Invasion" & alt == TRUE & outcome == "Military") {
    treat.name <- "Support for\n Use of Force"
    xlab <- "Indirect Effect of Invasion"
  }
  
  if(treatment.value == "Invasion" & alt == TRUE & outcome == "Taiwan") {
    treat.name <- "Support for Use of\n Force Against Taiwan"
    xlab <- "Indirect Effect of Invasion"
  }
  
  if(treatment.value == "Sanction") {
    treat.name <- "B. Economic Measures"
  }
  
  if(treatment.value == "Military") {
    treat.name <- "C. Military Measures"
  }
  
  if(treatment.value == "LackMilitary") {
    treat.name <- "D. Lack of Military Measures"
  }
  
  res2_1 <- out1$bin.result$results$indirect.effect$Group["est", ]
  c <- c(NA, NA, NA)
  names(c) <- c("y1.Legality", "y1.US_Threat", "y1.Peaceful_Image")
  res2_1 <- c(res2_1, c)
  res2_se_1 <- out1$bin.result$results$indirect.effect$Group["sd", ]
  res2_se_1 <- c(res2_se_1, c)
  
  
  if(treatment.value != "LackMilitary") {
    res2_2 <- out2$bin.result$results$indirect.effect$Group["est", ]
    res2_se_2 <- out2$bin.result$results$indirect.effect$Group["sd", ]
  }
  
  m2_1 <- data.frame(cbind(res2_1, res2_se_1))
  order <- order(m2_1$res2_1, decreasing = T)
  m2_1 <- m2_1[order, ]
  
  if(treatment.value == "LackMilitary") {
    res2_2 <- rep(NA, 11)
    res2_se_2 <- rep(NA, 11)
    names(res2_2) <- names(res2_se_2) <- names(res2_1)
    
    m2_2 <- data.frame(cbind(res2_2, res2_se_2))
    m2_2 <- m2_2[order, ]
  } else {
    m2_2 <- data.frame(cbind(res2_2, res2_se_2))
    m2_2 <- m2_2[order, ]
  }
  
  
  m2_1$term <- paste0("t", 1:nrow(m2_1))
  m2_2$term <- paste0("t", 1:nrow(m2_2))
  
  colnames(m2_1) <- c("estimate", "std.error", "term")
  rownames(m2_1) <- gsub("y1.", "", rownames(m2_1))
  
  colnames(m2_2) <- c("estimate", "std.error", "term")
  rownames(m2_2) <- gsub("y1.", "", rownames(m2_2))
  
  m2_1$model <- 0
  m2_2$model <- 1
  
  m2 <- rbind(m2_2, m2_1)
  
  if(outcome == "Military" & legend == 1) {
    g1 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "All Mediators", 
                             t2 = rownames(m2)[2], 
                             t3 = rownames(m2)[3], 
                             t4 = rownames(m2)[4],
                             t5 = rownames(m2)[5], 
                             t6 = rownames(m2)[6], 
                             t7 = rownames(m2)[7], 
                             t8 = rownames(m2)[8],
                             t9 = rownames(m2)[9], 
                             t10 = rownames(m2)[10], 
                             t11 = rownames(m2)[11]                             
        )) +
        theme_classic() + xlab(xlab) + ylab("Mediator") + xlim(-0.15, 0.30)+
        ggtitle(paste0("", treat.name)) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 18),
              axis.text.y  = element_text(size = 18),
              text = element_text(size = 18),
              legend.text = element_text(size = 18),
              legend.position = c(1.0, 0.22),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))      } } 
  if(outcome == "Military" & legend == 0) {
    g1 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "All Mediators", 
                             t2 = rownames(m2)[2], 
                             t3 = rownames(m2)[3], 
                             t4 = rownames(m2)[4],
                             t5 = rownames(m2)[5], 
                             t6 = rownames(m2)[6], 
                             t7 = rownames(m2)[7], 
                             t8 = rownames(m2)[8],
                             t9 = rownames(m2)[9], 
                             t10 = rownames(m2)[10], 
                             t11 = rownames(m2)[11]                             
        )) +
        theme_classic() + xlab(xlab) + ylab("Mediator") + xlim(-0.15, 0.30)+
        ggtitle(paste0("", treat.name)) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 18),
              axis.text.y  = element_text(size = 18),
              text = element_text(size = 18), legend.position = "none") + 
        scale_colour_grey(start = 0.1, end = .7)      } } 
  
  if(outcome != "Military" & legend == 1) {
    g1 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "All Mediators", 
                             t2 = rownames(m2)[2], 
                             t3 = rownames(m2)[3], 
                             t4 = rownames(m2)[4],
                             t5 = rownames(m2)[5], 
                             t6 = rownames(m2)[6], 
                             t7 = rownames(m2)[7], 
                             t8 = rownames(m2)[8],
                             t9 = rownames(m2)[9], 
                             t10 = rownames(m2)[10], 
                             t11 = rownames(m2)[11]
        )) +
        theme_classic() + xlab(xlab) + ylab("Mediator") + xlim(-0.15, 0.30)+
        ggtitle(paste0("", treat.name)) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 18),
              axis.text.y  = element_text(size = 18),
              text = element_text(size = 18),
              legend.text = element_text(size = 18),
              legend.position = c(1.0, 0.22),
              legend.justification = c(1, 1),
              legend.background = element_rect(colour="white"),
              legend.title.align = .5) + 
        scale_colour_grey(start = 0.1, end = .7,
                          name = "",
                          breaks = c("1", "0"),
                          labels = c("Experiment 2", "Experiment 1"))       }                          
    
  }
  if(outcome != "Military" & legend == 0) {
    g1 <- {dwplot(m2, 
                  dot_args = list(size = 2.5), whisker_args = list(size = 0.75),  
                  vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) %>% 
        relabel_predictors(c(t1 = "All Mediators", 
                             t2 = rownames(m2)[2], 
                             t3 = rownames(m2)[3], 
                             t4 = rownames(m2)[4],
                             t5 = rownames(m2)[5], 
                             t6 = rownames(m2)[6], 
                             t7 = rownames(m2)[7], 
                             t8 = rownames(m2)[8],
                             t9 = rownames(m2)[9], 
                             t10 = rownames(m2)[10], 
                             t11 = rownames(m2)[11]
        )) +
        theme_classic() + xlab(xlab) + ylab("Mediator") + xlim(-0.15, 0.30)+
        ggtitle(paste0("", treat.name)) +
        theme(plot.title = element_text(face="bold"), 
              axis.text.x  = element_text(size = 18),
              axis.text.y  = element_text(size = 18),
              text = element_text(size = 18), legend.position = "none") + 
        scale_colour_grey(start = 0.1, end = .7)      }                          
    
  }
  
  g1
}

## ----------------------
## Mediation Imai et al
## ----------------------

plots.mediation <- function(data, treatment, control, outcome, exp, filename) {
  subset <- data[which(data$RussianGroup %in% c(control, treatment) & is.na(data[, c(outcome)]) == F),]
  Dat_Mediate <- na.omit(subset)
  
  model.m.econ <- lm(Econ_Cost ~ Group, data = Dat_Mediate)
  model.y <- lm(paste(outcome, "~ Group + Econ_Cost"), data = Dat_Mediate)
  out.econ <- mediate(model.m.econ, model.y, treat = "Group",
                      mediator = "Econ_Cost")
  
  model.m.military <- lm(Military_Cost ~ Group, data = Dat_Mediate)
  model.y <- lm(paste(outcome, "~ Group + Military_Cost"), data = Dat_Mediate)
  out.military <- mediate(model.m.military, model.y, treat = "Group",
                          mediator = "Military_Cost")
  
  model.m.moral <- lm(Morality ~ Group, data = Dat_Mediate)
  model.y <- lm(paste(outcome, "~ Group + Morality"), data = Dat_Mediate)
  out.moral <- mediate(model.m.moral, model.y, treat = "Group",
                       mediator = "Morality")
  
  model.m.us <- lm(US ~ Group, data = Dat_Mediate)
  model.y <- lm(paste(outcome, "~ Group + US"), data = Dat_Mediate)
  out.us <- mediate(model.m.us, model.y, treat = "Group",
                    mediator = "US")
  
  model.m.peace <- lm(Peace ~ Group, data = Dat_Mediate)
  model.y <- lm(paste(outcome, "~ Group + Peace"), data = Dat_Mediate)
  out.peace <- mediate(model.m.peace, model.y, treat = "Group",
                       mediator = "Peace")
  
  model.m.success <- lm(Success ~ Group, data = Dat_Mediate)
  model.y <- lm(paste(outcome, "~ Group + Success"), data = Dat_Mediate)
  out.success <- mediate(model.m.success, model.y, treat = "Group",
                         mediator = "Success")
  
  model.m.threat <- lm(Threat~ Group, data = Dat_Mediate)
  model.y <- lm(paste(outcome, "~ Group + Threat"), data = Dat_Mediate)
  out.threat <- mediate(model.m.threat, model.y, treat = "Group",
                        mediator = "Threat")
  
  if(exp == 2) {
    model.m.legal <- lm(Legality ~ Group, data = Dat_Mediate)
    model.y <- lm(paste(outcome, "~ Group + Legality"), data = Dat_Mediate)
    out.legal <- mediate(model.m.legal, model.y, treat = "Group",
                         mediator = "Legality")
    
    model.m.usthreat <- lm(US_Threat ~ Group, data = Dat_Mediate)
    model.y <- lm(paste(outcome, "~ Group + US_Threat"), data = Dat_Mediate)
    out.usthreat <- mediate(model.m.usthreat, model.y, treat = "Group",
                            mediator = "US_Threat")
    
    model.m.peaima <- lm(Peaceful_Image ~ Group, data = Dat_Mediate)
    model.y <- lm(paste(outcome, "~ Group + Peaceful_Image"), data = Dat_Mediate)
    out.peaima <- mediate(model.m.peaima, model.y, treat = "Group",
                          mediator = "Peaceful_Image")
    
  }
  
  if(exp == 1) {
    pdf(file = paste0("./Figures_Appendix/", filename,".pdf"), width = 8, height = 6)
    par(mfrow = c(2, 4))
    plot(out.econ, xlim = c(-0.03, 0.35), main = "Economic Costs")
    plot(out.military, xlim = c(-0.03, 0.35), main = "Military Costs")
    plot(out.moral, xlim = c(-0.03, 0.35), main = "Morality")
    plot(out.us, xlim = c(-0.03, 0.35), main = "US Influence")
    plot(out.peace, xlim = c(-0.03, 0.35), main = "Peaceful Resolution")
    plot(out.success, xlim = c(-0.03, 0.35), main = "Likelihood of Success")
    plot(out.threat, xlim = c(-0.03, 0.35), main = "Threat")
    if(outcome == "Military"){
      mtext("The Effect of the Invasion Treatment on the Support for the Use of Force", side = 1, line = -24, outer = TRUE)
    }
    if(outcome == "Taiwan"){
      mtext("The Effect of the Invasion Treatment on the Support for the Use of Force Against Taiwan", side = 1, line = -24, outer = TRUE)
    }
    dev.off()
  } else {
    pdf(file = paste0("./Figures_Appendix/", filename, ".pdf"), width = 8, height = 9)
    par(mfrow = c(3, 4))
    plot(out.econ, xlim = c(-0.03, 0.35), main = "Economic Costs")
    plot(out.military, xlim = c(-0.03, 0.35), main = "Military Costs")
    plot(out.moral, xlim = c(-0.03, 0.35), main = "Morality")
    plot(out.us, xlim = c(-0.03, 0.35), main = "US Influence")
    plot(out.peace, xlim = c(-0.03, 0.35), main = "Peaceful Resolution")
    plot(out.success, xlim = c(-0.03, 0.35), main = "Likelihood of Success")
    plot(out.threat, xlim = c(-0.03, 0.35), main = "Threat")
    plot(out.legal, xlim = c(-0.03, 0.35), main = "Legality")
    plot(out.peaima, xlim = c(-0.03, 0.35), main = "Peaceful Image")
    plot(out.usthreat, xlim = c(-0.03, 0.35), main = "US Threat")
    
    if(outcome == "Military"){
      mtext("The Effect of the Invasion Treatment on the Support for the Use of Force", side = 1, line = -2, outer = TRUE)
    }
    if(outcome == "Taiwan"){
      mtext("The Effect of the Invasion Treatment on the Support for the Use of Force Against Taiwan", side = 1, line = -2, outer = TRUE)
    }
    dev.off()
  }
}






